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ABSTRACT 

in ■ 

I We present results from fully general relativistic three-dimensional numerical studies of thick- 

' disk accretion onto a rapidly-rotating (Kerr) black hole with a spin axis that is tilted (not aligned) 

^ , with the angular momentum vector of the disk. We initialize the problem with the solution for 

I an aligned, constant angular momentum, accreting thick disk, which is then allowed to respond 

'"^ ' to the Lense-Thirring precession of the tilted black hole. The precession causes the disk to warp, 

, beginning at the inner edge and moving out on roughly the Lense-Thirring precession timescale. 

' The propagation of the warp stops at a radius in the disk at which other dynamical timescales, 

^ ■ primarily the azimuthal sound-crossing time, become shorter than the precession time. At this 

, point, the warp effectively freezes into the disk and the evolution becomes quasi-static, except in 

' cases where the sound-crossing time in the bulk of the disk is shorter than the local precession 

■ timescale. We see evidence that such disks undergo near solid-body precession after the initial 

I warping has frozen in. Simultaneous to the warping of the disk, there is also a tendency for the 

midplane to align with the symmetry plane of the black hole due to the preferential accretion of 
the most tilted disk gas. This alignment is not as pronounced, however, as it would be if more 
efficient angular momentum transport (e.g. from viscosity or magneto-rotational instability) were 
considered. 

o ■ 

H ■ 

I Subject headings: accretion, accretion disks — black hole physics — hydrodynamics — methods: 

' numerical — relativity 

> : 

' 1. Introduction 

P.: 

Quasars, active galactic nuclei (AGN), X-ray binaries, core-collapse supernovas, and gamma-ray bursts 
(GRBs) are among the most luminous objects in the universe and are therefore of considerable observational 
interest. They are also all likely to contain a central, strongly gravitating, compact object and be powered by 
some form of accretion. They are therefore also of considerable theoretical interest. However, the nonlinear, 
time-dependent, multi-dimensional nature of the physical laws governing these objects makes a study of 
them from first principles difficult. While such studies can provide critical insight, there is strong motivation 
to apply a numerical approach. 

Many numerical simulations have been carried out over the past three decades studying the time- 
dependent evolution of accretion flows around black holes in the hydrodynamic (e.g. Wilson 1972; Hawley 
et al. 1984; Hawley 1991) and MHD (Koide et al. 1999; Gammie et al. 2003; De Villiers & Hawley 2003) 
regimes. However, due to the complex nature of this problem, a variety of simplifying assumptions are 
normally implemented. A common simplification in much of the early work was to assume a non-rotating 
(Schwarzschild) black hole, often approximated using a pseudo-Newtonian (Paczynsky & Wiita 1980) poten- 
tial. A notable exception was the work of Wilson (1972), which considered the spherical infall of material with 
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a non-zero specific angular momentum toward a Kerr black hole using the full metric, although restricted 
to two spatial dimensions. Although the assumption of a non-rotating black hole considerably simplifies 
the problem while still allowing useful insight, this does not appear to be well represented among observed 
black-hole systems (e.g. Elvis et al. 2002; Yu & Tremaine 2002), nor does it agree with studies of black hole 
spin evolution (Gammie et al. 2004; Volonteri et al. 2004) which suggest that most black holes are likely to 
be spinning rapidly. 

The spin of a black hole (often denoted by the specific angular momentum a = J/M) can have important 
consequences on the accretion rate and the geometry of the accretion flow. One obvious example is the 
dependence of Vms, the radius of the marginally stable circular orbit, on the spin parameter [r^s = Qro 
for a = (non-rotating black hole); rms = tq for a = +M (maximal progradc rotation); rms = Src for 
a = —M (maximal retrograde rotation), where re = GM/c^ is the gravitational radius (Bardeen et al. 1972)]. 
Another example is the reversal of flow direction in the vicinity of a rotating black hole for a retrograde 
accretion disk (see Figure 6 of De Villiers & Hawley 2002). 

Another common simpliflcation in black hole accretion models has been the assumption that the mid- 
plane of the accretion flow is aligned with the symmetry plane of the black hole. This has historically been 
necessitated in numerical work by limitations in computing resources which have often limited simulations 
to two-dimensional axisymmetry. The justification for this treatment often relies on the Bardeen-Petterson 
eff'ect (Bardeen & Petterson 1975) through which a spinning black hole may align the inner region of the 
accretion flow with the symmetry plane of the black hole. However, now that computational resources are 
available to perform three-dimensional simulations of black-hole accretion, there are many compelling rea- 
sons to revisit the question of black- hole - disk alignment. First, the formation avenues for many black- hole - 
disk systems favor, or at least allow for, an initially tilted configuration (see Fragile et al. (2001) for a review 
of the arguments). Second, although the Bardeen-Petterson effect may align the inner disk, this alignment 
may not extend nearly as far out as was originally expected [15 — SOrc (Nelson & Papaloizou 2000) rather 
than 100 — lOOOrc] . Many numerical simulations extend to larger radii than this and, therefore, may not 
be justified in assuming that the midplane of the disk is aligned with the symmetry plane of the black hole. 
Finally, although there is a residual torque between the tilted outer disk and black hole which will eventually 
completely align the system (Rees 1978; Scheuer & Feiler 1996; Natarajan & Pringle 1998), this happens 
over such a long timescale that these systems remain tilted for most of their observational lifetimes. 

If an accretion disk is misaligned or tilted, it will be subject to differential Lense-Thirring precession. 
This is ultimately what gives rise to the Bardeen-Petterson effect. For an ideal test particle in a slightly tilted 
orbit at a radius r around a black hole of mass M and specific angular momentum a, gravitomagnetic torques 
cause the orbit to precess at a frequency JIlt = 2aM/r^. Close to the black hole, this is comparable to the 
orbital frequency [^Kep = (M/r^)^/^]. However, because of its strong radial dependence, Lense-Thirring 
precession becomes much weaker far from the hole. Nevertheless, since it has a cumulative effect, it can build 
up appreciably over a sufficiently large number of orbits. In an accretion disk, this cumulative build-up is 
limited by dynamic responses within the disk. For a steady-state disk, Lense-Thirring precession is expected 
to be important out to a unique, nearly constant transition radius (Bardeen & Petterson 1975; Kumar & 
Pringle 1985). 

A comprehensive understanding of tilted accretion disks can be aided by the insights of direct nu- 
merical simulations. An important first step was taken by Nelson & Papaloizou (2000), who utilized a 
three-dimensional Newtonian SPH code to study tilted thin accretion disks subject to a post-Newtonian 
gravitomagnetic torque. In the current work we present results from three-dimensional, fully relativistic 
simulations of tilted black-hole accretion flows. A general relativistic treatment ensures that all important 
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relativistic features, such as the cusp in the potential and the Einstein and Lense-Thirring precessions of 
the orbits, as well as any higher order couplings of these features, are treated properly. In this first study, 
we consider relatively simple, inviscid, thick accretion disks. These disks are well studied in the literature 
and provide a convenient starting point for our work. We should mention, however, that effective angular 
momentum transport (such as through viscous coupling as considered in Nelson & Papaloizou (2000)), would 
likely have a considerable impact on our results. 

Unless otherwise stated, we use geomctrizcd units (G = c — 1) throughout this paper. Units of length 
are parameterized in terms of the gravitational radius of the black hole, re = GM/c^. We also adopt the 
standard convention in which Greek (Latin) indices refer to 4(3)-space components. 



Geometrically thick disks (often referred to as black hole tori or toroidal neutron stars) have been 
proposed to form through a number of scenarios, including the collapse of supermassive neutron stars (Vietri 
& Stella 1998), the iron-core collapse of a massive stars (MacFadyen & Woosley 1999), and the mergers of 
black hole - neutron star binaries (Lee & Kluzniak 1999). They may also exist in the inner regions of more 
standard accretion flows. These disks are convenient to study numerically as they represent an equilibrium 
solution in the limit that self-gravity, radiation, and viscosity are negligible. However, they are subject to 
numerous instabilities: 1) the so-called "runaway instability" in tori with non-negligible mass and accretion 
(Font & Daignc 2002); 2) the magneto-rotational instability (MRI, Balbus & Hawley 1991); and 3) the 
Papaloizou-Pringle instability (PPI, Papaloizou & Pringle 1984). In this work, we assume that the mass 
of the torus is negligible compared to the mass of the black hole, thus allowing us to neglect the runaway 
instability. This also allows us to ignore the self-gravity of the torus gas and treat the background metric as 
fixed. Wc also ignore magnetic fields in the current work, so the MRI is not applicable. Below we discuss 
how the PPI affects our work in more detail. 

The analytic solution for the structure of a thick disk has been worked out for both non-rotating and 
rotating black holes in the limit of axisymmetry (Fishbone & Moncrief 1976; Abramowicz et al. 1978). 
However, no such solution has been worked out for more general tilted flows. Furthermore, it is not clear 

that a steady-state solution would result. In light of the violent formation scenarios proposed above, however, 
it seems likely a newly formed torus will be highly disturbed and possibly have an orbital plane that is tilted 
with respect to the symmetry plane of the black hole. This, in fact, is part of the motivation for this 
numerical study. Nevertheless, we are left with a decision about how to initialize the problem. We choose 
to initialize all our disks with the analytic solution for an axisymmetric thick disk, which is subsequently 
allowed to respond to the Lense-Thirring precession of a tilted black hole. 

For untilted accreting thick disks, the gas flows in an effective (gravitational plus centrifugal) potential. 
Since axisymmetric thick disks have been described in great detail elsewhere (e.g. Kozlowski et al. 1978), 
we present only a brief discussion here. Surfaces of constant pressure in the torus are determined from the 
relativistic analog of the Newtonian effective potential 



where P is the fluid pressure, p is the density, h is the relativistic enthalpy, and is the potential at 
the boundary of the torus (where P = 0). For constant angular momentum I, the form of the potential 
reduces to $ = ln(— Ut), where Ut is the speciflc binding energy. Provided I > Ijns-, where l^s is the angular 
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momentum of the marginally stable Keplerian orbit, the potential $(r, 6') will have a saddle point ^cusp 
at r = Tcusp, d = it/ 2. We can define the parameter A$ = $j„ — $ciisp as the potential barrier (energy 
gap) at the inner edge of the torus. If A$ < 0, the torus lies entirely within its Roche lobe. Such a 
configuration is marginally stable with respect to local axisymmetric perturbations and unstable to low- 
order non-axisymmctric PPI modes (Papaloizou & Pringlc 1984). If A$ > 0, the torus overflows its Roche 
lobe and accretion occurs through pressure-gradient forces across the cusp. This accretion occurs without 
dissipation of angular momentum; thus viscosity is not required for accretion. This accretion generally 
suppresses the growth of instabilities in the torus (Hawley 1991). All of the simulations presented in this 
work begin with A$ > 0. 



3. Hydrodynamic Equations in Kerr Geometry 

In this work, we solve the system of equations for relativistic hydrodynamics in flux- conserving form 
(Anninos & Fragile 2003) 



dD diDV 



dt dx^ 



0, (2) 



dE d(EV') dW d(WV') 



dt dx'^ dt dx 



0, (4) 



at dx' 2S° dxi ^ ""dxi 

where g is the determinant of the 4-mctric, W = ^J—g%^ is the relativistic boost factor, D = W p is the 
generalized fluid density, P — {Y — 1)E/W is the fluid pressure, = u'^/u* is the transport velocity. 
Si = Wphui is the covariant momentum density, and E = We = Wpe is the generalized internal energy 
density. 

These equations are evolved in a "tilted" Kerr-Schild polar coordinate system {t, r, i?, ip). This coordinate 
system is related to the usual (untilted) Kerr-Schild coordinates {t, r, 6, (j)) through a simple rotation about 
the y-axis by an angle /3o, such that 

sin d cos ip \ / cos /3o — sin Po \ / sin 6 cos c6 \ 

sini?sin(p = 1 sin^sin^ 1 . (5) 

cosi? / \ sin/3o cos/3o / \ cos^ / 

The Kerr line element in the tilted {t, r, (p) coordinate system is 

, o / 2Mr\ , 2 4Mr , , 4Marsin2 6'„ , 
ds^ = - 1 5- dt^ + -^dtdr ^ PidMi? 

V r y Q e 

^2dMv + (^1 + j dr^ 

(2Mr\ 
1 - — 2- j Pidrdi? 

-2asin2 q{\- P-^drd^p 

+ {q^TI + API) di?2 + 2 (^^TiTs + ^PiPa) dM^p 

+ (^2T| + AP|)d<^^ (6) 
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(7) 

where sm9, cos9, sin0, and cost/) arc given by equation 5. For a Schwarzschild (a = 0) or an aUgned Kerr 
black hole (/3o = nn, where n = 0, 1, 2,...), the {t, r, i?, ip) and (t, r, 0, (^) frames are equivalent. 

The computational advantages of the "horizon-adapted" Kerr-Schild form of the Kerr metric are de- 
scribed in Papadopoulos & Font (1998) and Font et al. (1998). The primary advantage is that, unlike 
Boyer-Lindquist coordinates, there are no singularities in the metric terms at the event horizon. This is 
particularly important for numerical calculations as it allows one to place the grid boundaries inside the 
horizon, thus ensuring that they are causally disconnected from the rest of the flow. 



4. Numerical Methods 

Our main purpose in this work is to examine the structure of tilted accreting thick disks orbiting around 

a Kerr black hole. In practice, wc actually initialize the disks aligned with the computational grid and tilt 
the black hole through a transformation of the metric. We have tested problems with the disk tilted relative 
to the grid and found similar results. We carry out these simulations using Cosmos, a massively parallel, 
multidimensional, radiation-chemo-magnetohydrodynamic code for both Newtonian and relativistic flows. 
The relativistic capabilities and tests of Cosmos are discussed in Anninos & Fragile (2003). The present 
work utilizes the zone-centered artificial viscosity hydrodynamics package in Cosmos. 

The simulations are initialized with the analytic solution for an axisymmetric torus with constant 
specific angular momentum I around a rotating black hole with spin a. For the initialization, we assume an 
isentropic equation of state P = Kp^, although during the evolution the adiabatic form P = (F — 1)E/W is 
used to recover the pressure when solving equations (2) - (4). For all but one model, wc set F = 4/3 and 
K = 1.06 X 10^^ (in cgs units). The surface of the torus is defined by the effective potential at the boundary 
of the torus i>,„ , although in practice the torus is cut off where the fluid density and pressure match the 
background values. All of our runs are initialized with a positive energy gap A$ = — ^cusp > 0, where 
^cusp is the saddle point in the potential at r = rcusp, & = 7r/2. With this choice, a small fraction of the 
torus mass is initially out of equilibrium. 

The simulations are carried out on a grid extending from 0.98rBH < r < r^ax^ < d < tt, and 
< (fi < 2-jT, where vbh is the radius of the black hole horizon. We choose r^ax such that the disk is initially 
contained entirely within the computational domain. In order to increase the resolution in the inner region of 
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the disk, we replace the radial coordinate r with a logarithmic coordinate 77 = l + \n{r/rBH)- In terms of this 
logarithmic coordinate, the grid extends from 0.98 < 77 < ijmax, with a typical radial grid spacing of 0.2ra 
near the horizon and Gra near the outer boundary. Zones are evenly spaced in both angular dimensions. The 
grid is normally resolved with 96 x 64 x 96 zones along the dimensions r], ip. This resolution is sufficient 
to resolve key disk properties such as the PPL To illustrate this point and demonstrate the ability of our 
code to reproduce other published work, we plot in Figure 1 the growth of the m = 1 and m = 2 PPI modes 
for a simulation comparable to Model A3p of De Villiers & Hawley (2002) using a 64 x 64 x 64 grid. Our 
results are very similar to those of de Villiers & Hawley as seen by comparing to their Figure 7. Additional 
verification of our results is provided by comparing simulations at different grid resolutions in order to study 
the sclf-convcrgcncc of our solutions. 

In the "background" regions not determined by the initial torus solution, we initialize the gas following 
the spherical Bondi accretion solution. We fix the parameters of this solution such that the rest mass 
present in the background is negligible (< 10^^) compared to the mass in the torus. Similarly, the mass 
accretion from the backgroimd is generally small (< 10^'^) compared to the mass accretion from the disk. 
The outer radial boundary is held fixed with the analytic Bondi solution for all evolved fields. The inner 
radial boundary uses simple flat (zero gradient) boundary conditions. Data are shared appropriately across 
angular boundaries. 

After the torus and background are initialized, the black hole is tilted by an angle /3o through a trans- 
formation of the metric. Each model is evolved to t > 10To,.b, where Tort is the orbital period measured at 
the pressure (and density) maximum of the axisymmetric torus (denoted reenter)- Most of the disks achieve 
a quasi-steady state [i.e. variations in measured disk properties, such as the twist and tilt, are small (< 10%) 
over dynamical timescales] by the end of the simulations. The principle exceptions are disks that undergo 
near solid-body precession and disks that accrete a large percentage of their mass (> 50%) onto the black 
hole. 

To describe our results more thoroughly, we track a number of quantities throughout the simulations, 
including the mass accretion rate, the PPI mode growth for the azimuthal wavenumbers m = 1 and 2, and 
the precession angle (twist) and angle of inclination (tilt) of the disk. The mass accretion rate is defined as 



The accretion rates are converted to units with a polytropic constant k = 1.5 x 10^° (cgs units) and mass 
M = IMq, consistent with Igumenshchcv & Bcloborodov (1997). This value is then normalized by the 
Eddington accretion rate MEdd = ^Edd/c^ ~ 1-4 x 1O^^(M/M0) g s~^. The conversion of M is accomplished 
through the relation 



Using this common scaling allows for easy comparison of our results. Nevertheless, we caution the reader 
against attaching any particular physical significance to these values since the mass accretion rates and 
densities in these simulations are essentially arbitrary. 

The m = 1 and 2 Fourier modes of the PPI are extracted by computing azimuthal density averages of 
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Fig. 1. — PPI mode growth for comparison with Model A3p of De Villiers & Hawley (2002). 
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the following form (De Villiers & Hawley 2002) 

^e[km{r)]= / p{r,n/2,ip)cos{m(p)dip , (10) 
Jo 

lva[km{r)\= I p{r,TT/2,ip)sm{rrvp)d(p . (11) 
Jo 

The power in mode m is then 

fm = ^ T"""^ ln({Re[fc„(r)]}2 + {Im[fc„(r)]}2)dr , (12) 

where Vmin and Vmax are the approximate inner and outer edges of the disk. 
We define the precession angle (twist) as (Nelson & Papaloizou 2000) 



j{r) = arccos 
where 



■y 



(13) 



3bh = {-aM sin l3ox,0,aM cos Poz) (14) 
is the angular momentum vector of the black hole and 

^Diskir) = [{JDisk)lX, {JDisk)2y, {Joiskhz] (15) 

is the angular momentum vector of the disk in an asymptotically flat space, where 
and 

LI"" = J (a;T"° - x^T^"^) d^a;, (17) 

with S" = j T^^d?x, T^^^ = phu'^u" + Pg'^'', and Po is the initial tilt between the angular momentum vectors 

of the black hole and the disk. The equations for L^'^ and S'^ arc integrated over concentric radial shells of 
the grid. The unit vector y points along the axis about which the black hole is initially tilted and z points 
along the initial symmetry axis of the disk. Thus, from Equation (13), -f{r) = throughout the disk at 
t = 0. In order to capture twists larger than 180°, we also track the projection of 3bh x Joiskir) onto x, 

allowing us to break the degeneracy in arccos. Wc use the condition 7 = 1 to define a transition radius rx 
between an inner twisted disk (precession greater than 1 radian) and the outer undisturbed disk. 

The angle of inclination (tilt) in the disk is defined as (Nelson & Papaloizou 2000) 



l3{r) = arccos 



l^BHllJoiskir) 



(18) 



The models are set up such that, initially, (3{r) = /3o throughout the disk. For tilted disks (/3o 7^ 0) such as 
the ones considered in this work, Lense-Thirring precession proceeds most rapidly closest to the black hole 
and gradually affects the disk further out as the simulations progress. 

We expect differential Lense-Thirring precession to dominate whenever the precession timescale tlt = 
Q^T shorter than local dynamical timescales in the disk (Bardeen & Petterson 1975; Kumar & Pringle 
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1985). Wc consider two possible limiting timescales: the mass accretion timescale Tacc = (r — rBH)/V^ , 
where V is the average inflow velocity of the torus gas, and the azimuthal sound-crossing time defined as 
Tcs = Trr/cs, where Cg is the local sound speed in the density midplane. In this work we determine the 
accretion timescale numerically from V = M/X, where A(r) = J J OMdip. Initially the sound speed in the 
disk is given by = (F — 1)(1 — Ut/ut.in)- During the evolution, the local soimd speed is recovered from 
the fluid state through the relation = r(r — l)P/[(r — l)p + rP]. Since the sound speed is defined in the 
proper frame of the fluid, it is not strictly accurate to compare Tcs to quantities defined using the coordinate 
time (such as tlt)- However, the corrections are only of order Mr/{r'^ + a^) and should therefore be small 
(< 20%) for the radii of interest here. 

5. Results 

In this work we present results from eleven different models, varying the spin of the black hole, the 
angular momentum and tilt of the torus, the magnitude of the energy gap, the sound speed within the disk, 
and the size and resolution of the grid. The basic parameters for each model are listed in Table 1. Many 
of our parameters closely match those of Models II and III in Igumenshchev & Beloborodov (1997), thus 
providing us with an avenue to compare our results with other numerical simulations. However, we emphasize 
that there are a number of significant differences between their simulations and ours: all of the simulations 
reported here arc performed in three dimensions and all but one arc tilted; we use the Kerr-Schild form of 
the Kerr metric; we use a logarithmic radial coordinate for higher resolution near the black hole; and we 
choose a different inner surface potential $j„. In the following sections we describe the results of each of our 
models in more detail. 



5.1. Model t90 

This is the only model without an initial tilt between the symmetry plane of the black hole and the 
midplane of the torus. The black hole has a spin a/M = 0.9 and radius tbhItg = 1-44; the torus has 
a specific angular momentum l/M = 2.6088 and a density center at rcenter/ra = 3.45. The cusp in the 
potential is located at rcusp/rc = 1.76, 9 = 7r/2 with ^cusp = —0.04. We choose an inner surface potential 
^in = —0.01 giving an energy gap A$ = 0.03. The outer boundary is set at Tmax/rG = 80 and the simulation 
is evolved until t = 490M « lOrorb- This untilted simulation is included primarily to illustrate that our base 
model is not susceptible to the Papaloizou-Pringle instability, owing to the relatively large, positive energy 
gap at the cusp leading to a moderately high accretion rate that suppresses the growth of PPI (Hawley 
1991). In Figure 2a we show (a) the final distribution of torus gas density, (&) the time evolution of the mass 
accretion rate at the black hole horizon, and (c) the growth rates for the m = 1 and m = 2 PPI modes. 

The mass accretion rate (Figure 2a6) starts at a small value fixed by the background Bondi solution. 
Gas from the disk begins its infall from r > rcusp and so experiences a short delay before reaching the black 
hole horizon. After a short transition time (< 3t„ri,). M asymptotes toward a stationary value. However, 
the initial adjustment of the disk generates an outgoing wave within the disk that is reflected at the outer 
disk edge. The arrival of this wave back at the inner edge of the disk causes the "bump" in the accretion 
rate at late time. A similar behavior is seen in many of our models. Although this feature shows up vc;ry 
clearly in the accretion rate, it does not appear to have a significant effect on any of the other disk properties 
we track. This is consistent with our findings throughout this study that the tilt and twist of these disks 



Fig. 2a. — Model t90. (a) Linear density plot a.t t = 490Af sa lOrorb- For rendering purposes, the surface 
of the torus is set at psurf ~ O.lpmax and half the torus is cut away. To give a clearer view of the central 
region, the image is zoomed in to a small fraction of the total computational volume. The image is plotted 
in logarithmic coordinate space, where x = 77 sin cos y = rj sin's simp, and z = rycosi?. The torus is only 
initialized in the region r > Vcusp- The black hole (not pictured) lies at the center of the image. (6) Mass 
accretion rate M at the black hole horizon in units of the Eddington mass accretion rate MEdd (assuming 
Mbh ~ 1-A^0 and k — 1.5 x 10^" in cgs units), (c) PPI mode growth. 
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Fig. 2c.— 
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are not strongly responsive to the mass accretion rate. Furthermore, we note from our longer runs below 
(see Model t915L) that this feature does not recur at later times; after the wave reaches the inner edge of 
the disk, it passes through the cusp and black-hole horizon and is not reflected. It, therefore, has no further 
impact on the problem. Based upon this behavior, it appears that this feature is either not related to the 
axisymmetric oscillation modes studied in Zanotti ct al. (2003) or that the large value of A$ (energy gap) 
prevents the oscillation mode from being sustained. 

Although the PPI is seeded with small (< 1%) random fluctuations in the initial mass and internal 
energy densities, its growth is suppressed and the torus becomes essentially smooth and axisymmetric. This 
is consistent with the findings of De Villiers & Hawley (2002) that the PPI is suppressed whenever the 
accretion rate is high as in our models. In fact, this motivated our choice of parameters for these models, as 
we wanted to avoid the complication of PPI in this study. 

5.2. Models t915XL, t915L, & t915H 

These models begin from exactly the same initial conditions as Model t90 except that the black hole is 
tilted by an angle /3o — 15°. Thus, these disks are subject to Lense-Thirring precession. Model t915XL uses 
a resolution of 24 x 24 x 24 zones and Model t915L uses a resolution of 48 x 24 x 48; both are evolved until 
t = 1400M « 30rorb. Model t915H uses a resolution of 96 x 64 x 96 and is evolved untfl t = 490M « lOTorb- 
Figure 3a shows (a) the final distribution of torus gas density for Model t915H, (6) the cumulative precession 
7 and (c) tilt /3 of the disk as a function of the logarithmic radial coordinate rj, (d) the mass accretion rate 
M at the black hole horizon as a function of time, and (e) the time evolution of the transition radius tt 
defined as the radius where the twist equals 1 radian (7 = 1). It is clear from Figures 3aa, 3a&, and Sac 
that the disk is strongly twisted and warped at the end of the simulation. Nevertheless, the disk reaches a 
quasi-steady state before the simulations are stopped. This claim is supported by the following observations: 
1) in Figure 3a6, the twist profile changes only slightly from t = lOTorb to t = SOTorb', 2) in Figure 3ad, the 
mass accretion rate asymptotes toward a constant value; and 3) in Figure 3ae, the transition radius also 
asymptotes toward a constant value. 

Although the transition radius in Figure 3ae starts out closely following the expected growth of Lense- 
Thirring precession for ideal test particles (vt oc t^^^), it eventually approaches a limiting value of rx ^ Tr^. 
This is approximately coincident with the radius at which the local azimuthal sound-crossing time Tcs (dashed 
line in Figure 3a e) becomes shorter than the Lense-Thirring precession timescale t^t {thick solid line). This 
finding appears consistent with the expectations of Nelson & Papaloizou (2000) that warps in thick disks are 
diffused by wave modes. As the disk is twisted up by Lense-Thirring precession, the twist is diffused away 
by dynamical responses in the disk. We note that the accretion timescale Tqcc (dot-dashed line) in these 
disks is significantly longer than the Lense-Thirring or sound- crossing timescales except in the outer regions 
(r > ISrc)- As such, it does not appear to play a role in limiting the precession. Nevertheless, we point out 
that this tilted disk is subject to a higher accretion rate than the untilted model (compare Figures 3a.d and 
2a6). This is expected since more of the tilted torus begins the simulation out of equilibrium. 

A curious feature seen in Figure 3ac and evident in nearly all of our models is the tendency for these 
disks to align toward the symmetry plane of the black hole despite the lack of viscous angular momentum 
transport. This alignment seems to be facilitated by the preferential accretion of highly tilted disk material. 
Notice in Figure 3aa that two opposing accretion streams start from latitudes well above and below the 
black-hole symmetry plane. The start of the lower stream is visible near the bottom of the disk. The upper 
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Fig. 3a. — Model t915L & t915H. (a) Linear density profile (similar to Figure 2aa) for Model t915H at 
t = 490M « lOTorb- The black hole spin axis is tilted 15° to the left of the ^-axis. (6) Cumulative precession 
angle 7 and (c) tilt f3 in radians as functions of radius. Initially, 7 = and /3 — 15° « 0.26 throughout 
the disks, (d) Mass accretion rate in units of the Eddington mass accretion rate, (e) Transition radius 
(symbols) as a function of time along with the Lense-Thirring precession timescale t^t (solid), azimuthal 
sound crossing time Tcs (dashed), and accretion timescale Tacc (dot- dashed). The thin vertical line marks the 
radius where tlt = Tcs- 
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stream actually begins in the portion of the disk that is cut away, but this stream is nevertheless prominently 
visible cutting across the upper half of the disk. These streams begin near the density center of the disk and 
are 180° out of phase with one another. The material in these streams has a very large tilt relative to the 
black hole and has a larger tilt, on average, than the material slightly further out in the disk. This explains 
the prominent increase in tilt toward smaller radii seen near reenter in Figure Sac. Closer to the black hole, 
these two accretion streams spiral around until they cross near the symmetry plane of the black hole. As 
they cross the coUisional dissipation of angular momentum produces a small ring of gas orbiting very near 
the symmetry plane of the black hole. This ring has a very small tilt, on average, and is responsible for the 
small value of /? seen near the black hole horizon in Figure Sac. 

The preferential accretion into the black hole of material with a large average tilt removes much of the 
misaligned angular momentum from the disk, thus allowing the system to gradually align. For the low-mass 
disks considered here, it is primarily the disk orientation that changes. However, for sufficiently massive 
disks, such accretion could gradually change the orientation of the black hole spin axis. 

We note from all of the figures that there is good agreement between the runs at different resolutions. 
This is both reassuring and convenient as it allows us to use the course grid to follow the evolution of 
these models for longer periods at less computational expense. We also note that the results appear to be 
converging since the results of our intermediate resolution model (t915L) agree more closely with those of 
the high resolution model (t915H) than the low resolution one (t915XL). Nevertheless, we acknowledge that 
there may be features which are still not resolved even with our finest grid. 

5.3. Model t915H2 

The initialization of Model t915H2 differs from Model t915H above primarily in that Model t915H2 has 
a larger initial energy gap at the cusp (A$ = 0.04). The larger energy gap implies a higher mass accretion 
rate and a larger initial torus. The torus (and the grid) for Model t915H2 are about three times larger than 
in Model t915H {rmax/fG = 240). This model is resolved with the same high-resolution grid as Model t915H 
(96 X 64 X 96 zones), so the effective resolution in the radial direction is reduced by about a factor of three. 

Despite the higher accretion rate of this model (compare Figures 4ad and Sad), the measured precession 
7 (Figure 4a6) and tilt /3 (Figure 4ac) are very similar to those of Models t915L and t915H. The final density 
profile is also very similar (compare Figures 4aa and Saa). This suggests that the accretion rate is not a 
dominant factor controlling the evolution of these models. The size of the torus also does not appear to play 
a significant role. 

5.4. Models t930L & t930H 

These models also begin from exactly the same initial conditions as Model t90 except that the black 
hole is tilted by an angle /3a = S0°. These models allow us to gauge how our results depend on the initial tilt 
of the disk. Model t930L uses a resolution of 48 x 24 x 48 zones and is evolved until t — 1400M « SOrorb; 
Model t930H uses a resolution of 96 x 64 x 96 and is evolved until t = 460M w lOTorb- Again there is very 
good agreement between the two. 

These models have the highest accretion rates (Figure Sad) of any considered. Again, this is expected 
since they begin with the disk furthest from equilibrium. One consequence of the very high accretion rate is 
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Fig. 4a. — Model t915H2. (a) Linear density profile (similar to Figure 2aa) at t — 470A/ w IOtq,;,. The 
black hole spin axis is tilted 15° to the left of the ^-axis. (6) Cumulative precession angle 7 and (c) tilt /3 in 
radians as functions of radius. Initially, 7 = and /? = 15° « 0.26 throughout the disks, (d) Mass accretion 
rate in units of the Eddington mass accretion rate, (e) Transition radius (symbols) as a function of time 
along with the Lense-Thirring precession timescale tlt, azimuthal sound crossing time Tcs, and accretion 
timescale Tacc- The thin vertical line marks the radius where t^t — Tcs- 



- 21 - 



r (GM/c^) 
5 10 20 




Fig. 4b.— 



- 22 - 




- 23 - 




Fig. 4d.— 



-24- 




Fig. 4e. — 



- 25 - 



that the disk is largely depleted of gas by the end of the simulation, particularly for the longer evolution of 
Model t930L, which was carried out to t = 1400M « SOTorb- Nevertheless, we find very similar results for 
the twist and tilt of these f3o = 30° models compared with the previous /3o = 15° models (compare Figures 
5ac and Sac). We again see the preferential accretion of highly tilted gas through two opposing streams 
starting near reenter- Again these streams cross near the symmetry plane of the black hole, leaving a small, 
nearly aligned ring of gas close to the horizon. These results further confirm that the mass accretion rate 
plays little, if any, role in regulating the precession or tilt of these disks. 

As in the previous models, the transition radius tt (Figure 5ae) asymptotes toward a radius tt ^ Gtg 
consistent with the radius at which tlt ~ Tcs- The early migration of Vt is even similar to the previous 
models (compare Figures 5ae and 3ae). This leaves the larger accretion rate and tilt as practically the only 
differences between this and previous models. Thus, it appears our results are only weakly dependent on 
the initial tilt. 



5.5. Model t915r 

One of our general conclusions to this point is that differential Lense-Thirring precession in the disk is 
limited to radii for which tlt < Tcs- However, this conclusion is based on a series of models that all share the 
same precession frequency and roughly the same sound speed in the disk. It is useful, therefore, to consider 
a model where at least one of these timescales is changed. In Model t915r, we use the same black-hole spin 
{a/M = 0.9), disk parameters {l/M = 2.6088 and = -0.01), and resolution (48 x 24 x 48) as Model 
t915L, except we change the equation of state of the gas such that F = 5/3 and k = 1.24 x 10^ (in cgs units). 
Since the sound speed in these models scales as (F — 1)^/^, this change increases the sound speed in the disk 
by \/2 and decreases the azimuthal sound- crossing time by the same factor. 

All other things being equal, our first expectation is that the faster sound speed in this disk should make 

the transition radius asymptote at a smaller value than in the previous models. Instead we sec in Figure 6a6 
that the twist continues to grow with time throughout the disk. This is also apparent in Figure 6ae where 
the transition radius tt gradually moves outward in the disk throughout the entire simulation and does not 
asymptote toward a constant value. By the end of the simulation it reaches a value more than double the 
radius where t^t = Tcs- Although this result seems surprising in light of what we found in the previous 
models, when considered within the context of the results presented below, it is apparent that Model t915F 
is a transition model. In all of the previous models, the radius where tlt = Tcs lies outside of reenter', thus 
differential Lense-Thirring precession dominates the bulk of the disk material. In the models presented in 
the next two sections, the radius at which tlt = ^cs hes inside of reenter, those disks, therefore, have very 
efficient internal communication throughout the bulk of the disk and respond more uniformly to precession. 
For the current model, reenter and the radius at which tlt = Tcs differ by less than 25%; as a consequence, 
Tlt ~ Tcs for the bulk of the disk material. This balance prevents either differential precession or wave 
propagation from dominating the evolution. As a result, in Figure 6ae we find an intermediate behavior 
between the vertical asymptote seen in strong differential precession models (Figures 3a e, 4a e, and 5a e) and 
the nearly horizontal limit of solid body precession (Figure 8a e). 



- 26 - 




Fig. 5a.— Model t930L & t930H. (a) Linear density profile (similar to Figure 2aa) for Model t930H at 
t — 460Af « IQTorb- The black hole spin axis is tilted 30° to the left of the z-axis. (6) Cumulative precession 
angle 7 and (c) tilt [3 in radians as functions of radius. Initially, 7 = and (3 — 30° w 0.52 throughout the 
disks, (rf) Mass accretion rate in units of the Eddington mass accretion rate, (e) Transition radius (symbols) 
as a function of time along with the Lense-Thirring precession timescale tlt, azimuthal sound crossing time 
Tcs, and accretion timescale Tacc- The thin vertical line marks the radius where tlt — Tcs- 
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Fig. 6a. — Model t915r. (a) Linear density profile (similar to Figure 2aa) at t = 1400Af w SOrorb- The 
black hole spin axis is tilted 15° to the left of the z-axis. (&) Cumulative precession angle 7 and (c) tilt 
P in radians as functions of radius. Initially, 7 = and P = 15° « 0.26 throughout the disks, (d) Mass 
accretion rate scaled to allow direct comparison with our other models, (e) Transition radius (symbols) as a 
function of time along with the Lense-Thirring precession timescale t^t, azimuthal sound crossing time Tcs, 
and accretion timescale Tacc- The thin vertical line marks the radius where t^t = Tcs- 
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5.6. Model t915R 

Model t915R explores a retrograde torus model. The magnitude of the black hole spin (|o| = 0.9) is the 
same as in the previous models, but it spins in the opposite sense. Thus, the angular momentum of the black 
hole and torus are almost anti-parallel; they differ from anti-parallel by 15°. The retrograde orbit requires 
the torus to have a higher specific angular momentum to ensure / > Ims- In this case, l/M = 4.4751; the 
torus center is at Vcenter/fG = 12. The cusp in the potential is located at Vcusp/i'G = 6.67, 9 = 7r/2 with 
^cusp = —0.03. We choose an inner surface potential = —0.0125 giving an energy gap A4> = 0.0175. 
The outer boundary is set at rmax/rc = 80 and the simulation is evolved until t = 2000M « STorb- 

The most obvious thing we notice in Figure 7aa is that the accretion stream in this model takes the form 
of a relatively thin disk lying very close to the symmetry plane of the black hole. In this sense, the results 
appear very similar to the usual Bardeen-Petterson configuration, but for a very different physical reason as 
we shall see. The gas in the accretion stream still orbits in a retrograde sense just as the rest of the disk does. 
Therefore, the gas in the accretion stream has a tilt /? « tt (Figure 7ac) consistent with retrograde motion 
near the symmetry plane of the black hole. Closest to the black hole (ry < 1.6; r/rc < 2.62), the accreting 
gas actually appears to tilt away from the symmetry plane. This is due the inclusion of background gas in 
our calculations. Close to the black hole, the background gas is undergoing very rapid progradc precession 
near the symmetry plane of the black hole. It, therefore, provides a /3 « contribution near the black hole. 

It is important to note that, within the accretion stream at r < Vcusp, the accretion timescale Tqcc is 
actually shorter than both the sound-crossing time Tcs and the precession timescale tlt (Figure 7ae). We, 
therefore, don't expect precession to have a dramatic effect at these small radii in this model since the Lense- 
Thirring precession timescale is longer than the dynamical timescale. The role of Lense-Thirring precession 
inside of rcusp is further diminished by the fact that most of the gas inside this radius lies very close to the 
symmetry plane of the black hole making precession essentially meaningless. Notice, for instance, that the 
calculation of the precession 7, which relies on the cross-product of J bh and J Disk , which becomes unreliable 
whenever Jbh and J Disk are nearly parallel or anti-parallel. 

Outside of rcusp, the accretion timescale becomes significantly longer, yet the sound-crossing time quickly 
becomes shorter than the precession timescale. Therefore, true differential precession is only expected to be 
important over a very small region in this disk model. Importantly though, we see in Figure 7a6 that the 
precession proceeds in a retrograde fashion as expected (Lense-Thirring precession always proceeds in the 
direction of the black- hole spin). Precession also extends out to much larger radii (essentially throughout 
the disk) in this model than in previous models (see Figure 7a,b). We will return to this point in the next 
section. 



5.7. Model t515L & t515H 

Model t515 explores the dependence of our results on the magnitude of the black hole angular mo- 
mentum. For these models, the black hole has a spin a/M = 0.5 and radius tbhItg = 1-87; the torus 
has a specific angular momentum l/M = 3.385; the torus center is at rcenter/ra = 6.92. The cusp in the 
potential is located at Vcusp/i'G = 2.96, 6 = 7r/2 with ^cusp = —0.0136. We choose an inner surface potential 
$i„ = —0.0096 giving an energy gap A<I> = 0.004. The outer boundary is set at rmax/rc = 120. Model 
t515L uses a resolution of 48 x 24 x 48 zones and is evolved until t = 2400M « 2QTorb', Model t515II uses 
a resolution of 96 x 64 x 96 and is evolved until t = 1200M w lOTorb- Since the grid in these models has a 
larger maximum radius than most of the other models, the effective resolution is somewhat less. 



Fig. 7a. — Model t915R. (a) Linear density profile (similar to Figure 2aa) at t = 2000Af w Srorb- The 
black hole spin axis is tilted 15° to the right of the — z-axis (points down and to the right). (6) Cumulative 
precession angle 7 and (c) tilt (3 in radians as functions of radius. Initially, 7 = and (3 = 165° w 2.88 
throughout the disks, (d) Mass accretion rate in units of the Eddington mass accretion rate, (e) Transition 
radius (symbols) as a function of time along with the Lense-Thirring precession timescale tlt, azimuthal 
sound crossing time Tcs, and accretion timescale Tacc- Here the thin vertical line marks r = Vcusp- 
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There are two striking features to these models. First, in Figures 8aa and c, we see that the inner 
accretion disk appears to tilt away from the symmetry plane of the black hole. This is actually consistent 
with what we have seen in previous models in the sense that the black hole preferentially accretes the most 
tilted material. As in previous models, two opposing accretion streams begin from latitudes well above and 
below the black- hole symmetry plane near reenter- The accretion streams in this model actually begin at 
comparable latitudes to what is seen in Model t915 (compare Figures 8aa and 3aa). Again, this material 
carries much of the tilted angular momentum into the hole and allows the disk to gradually align. Since 
the streams have a greater tilt on average than the material further out, /3 increases toward smaller radii 
inside reenter as seen in Figure Sac. However, unlike our a/M = 0.9 models, these accretion streams do 
not cross before passing into the hole and we do not see the formation of an inner ring of low j3 material. 
Consequently, /? remains large from near reenter all the way in to the horizon. 

The second striking feature is apparent in Figure 8ae. Although the transition radius at first approaches 
an asymptotic value coincident with tlt ~ Tcs similar to our other models, it later appears to rapidly move 

outward with no sign of another limit being reached. This is actually evidence of near rigid-body precession 
in the disk. For pure rigid-body precession, the twist 7 in Figure 8a6 would be a horizontal line and the 
transition radius would lose its meaning since there would be no "transition" from an inner, twisted region 
to an outer, untwisted one. Model t515 doesn't quite undergo pure rigid-body precession. It actually starts 
o\it undergoing differential precession and then switches to near rigid-body precession. This change is most 
apparent in the time sequence of 7-profiles shown in Figure 8a6. This also explains the near constant 
precession angle for gas beyond r = reenter in Figures 7a 6. 

Generally, rigid-body precession is expected in disks whenever Tes is less than the applicable precession 

timescale (Larwood et al. 1996), in this case whenever Tes < tlt- This condition is actually met in all of our 
models at large enough radii; however, only in Models t515 and t915R is this condition met for r < reenter- 
Thus, only in these two models does the bulk of the gas in the disk meet the condition for rigid-body 
precession. 

6. Discussion and Summary 

We have performed the first fully rclativistic numerical studies of tilted accretion disks around Kerr 
black holes. Such disks are subject to Lense-Thirring precession, resulting in a torque that tends to twist 
and warp the disk. For the thick disks considered here, we find that the nature of this precession depends 
primarily on the somid speed in the disk. Whenever t^t < Tcs in the bulk of the disk, it undergoes differential 
precession out to a transition radius set by tlt ~ Tes ^ Tacc, which occurs at a relatively modest distance 
of r/rc ~ 7. In other disk models, such as viscous Keplerian disks, this transition may occur at somewhat 
larger radii (~ 15 — SOrc, Nelson & Papaloizou 2000). We find that the location of this transition radius 
does not depend strongly on the size, mass, mass accretion rate, or tilt of the disk. The initial tilt does, 
however, affect the overall mass accretion rate for our models, since higher tilt models begin further from 
equilibrium within the black hole potential. 

Whenever tlt > Tcs for the bulk of the disk, we find that it undergoes near rigid-body precession after a 
short initial period of differential precession. Such rigid-body precession of a thick disk could have important 
astrophysical consequences. For instance, this is the basis of at least one model attempting to explain the 
106 day variability observed from Sgr A* (Liu & Melia 2002), the 2.6 x IO^Mq black hole at the Galactic 
center. It could also be important for understanding the jet precession observed in Galactic microquasars 



Fig. 8a. — Model t515. (a) Linear density profile (similar to Figure 2aa) a,t t — 1200Af « lOTorb- The 
black hole spin axis is tilted 15° to the left of the z-axis. (6) Cumulative precession angle 7 and (c) tilt /3 in 
radians as functions of radius. Initially, 7 = and /3 = 15° « 0.26 throughout the disks. In panel 6, we have 
included data from t — 0.3, 1, 1.6, 3.2, 6.4, and Srorb {dotted lines) to give a more complete time history of 
the precession, (d) Mass accretion rate in units of the Eddington mass accretion rate, (e) Transition radius 
(symbols) as a function of time along with the Lense-Thirring precession timescale t^t, azimuthal sound 
crossing time Tcs, and accretion timescale r^cc- The thin vertical line marks the radius where tlt — Tcs- 
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such as SS 433 and Her X-1. 

One obvious shortcoming of this work is that we have ignored angular momentum transport mechanisms 
that are likely present in most astrophysical disks [e.g. viscosity or the magneto-rotational instability (Balbus 
& Hawley 1991)]. Angular momentum dissipation may allow the inner region of the tilted disk to align more 
quickly with the symmetry plane of the black hole than what we see in the inviscid disks considered here 
(Bardccn & Petterson 1975; Kumar & Pringle 1985; Schcuer & Feilcr 1996). The alignment may also be more 
efficient in viscous disks, resulting in a sharper transition (Nelson & Papaloizou 2000). A sharp transition 
from an inner aligned disk to a tilted outer disk may be useful in explaining a variety of phenomena observed 
in black-hole systems, such as quasi-periodic oscillations (Fragile et al. 2001) as well as certain spectral 
features. 

The authors would like to thank the Visit development team at Lawrence Livermore National Laboratory 
(http://www.llnl.gov/visit/), in particular Hank Childs, for visualization support. P.C.F. would like to 
acknowledge many useful discussions on this work, particularly with G. Mathews, J. Wilson, and C. Gammie. 
This work was performed under the auspices of the U.S. Department of Energy by University of California, 
Lawrence Livermore National Laboratory under Contract W-7405-Eng-48. 



Tabic' 1. 



Tilted Thick-Disk Models 



Grid 



Model 




a/M 


l/M 




A$ 


^0 


t90 


96 X 64 X 96 


0.9 


2.6088 


46.1 


0.03 





t915XL 


24 X 24 X 24 


0.9 


2.6088 


46.1 


0.03 


15° 


t915L 


48 X 24 X 48 


0.9 


2.6088 


46.1 


0.03 


15° 


t915H 


96 X 64 X 96 


0.9 


2.6088 


46.1 


0.03 


15° 


t915H2 


96 X 64 X 96 


0.9 


2.6088 


46.1 


0.04 


15° 


t930L 


48 X 24 X 48 


0.9 


2.6088 


46.1 


0.03 


30° 


t930H 


96 X 64 X 96 


0.9 


2.6088 


46.1 


0.03 


30° 


t915r'^ 


48 X 24 X 48 


0.9 


2.6088 


46.1 


0.03 


15° 


t915R 


96 X 64 X 96 


-0.9^^ 


4.4751 


255. 


0.0175 


15° 


t515L 


48 X 24 X 48 


0.5 


3.385 


118. 


0.004 


15° 


t515H 


96 X 64 X 96 


0.5 


3.385 


118. 


0.004 


15° 



"^Orbital period measured at the pressure center of the torus. 
•T = 5/3. 
"^Retrograde torus. 
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